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This paper is part of a program of investigation of the chiral transition in Nf = 2 QCD, started 
in Ref. [lj. Progress is reported on the understanding of some possible systematic errors. A direct 
test of first order scaling is presented. 



I. INTRODUCTION 

The order of the chiral transition of two-flavor QCD is 
not established yet, despite its great relevance to under- 
stand important aspects of color confinement and of the 
structure of the QCD phase diagram 0]. A natural order 
parameter in that limit is the chiral condensate, even if 
empirically one finds that also deconfinement happens at 
the same critical temperature T c , as found by looking at 
the susceptibility of the Polyakov loop or by direct stud- 
ies of order parameters constructed in the framework of 
specific models of color confinement 0, 0, H[ • 

A renormalization group analysis plus e-expansion can 
be made around m ~ 0, assuming that the relevant de- 
grees of freedom for the chiral transition are scalar and 
pseudoscalar fields [1,0, [1]. The result is that for Nf = 2, 
if the 17a (1) anomaly term is still effective at the transi- 
tion, i.e. if the rf mass does not vanish at T c , an infrared 
stable fixed point may exist in the universality class of 
the three dimensional 0(4) spin model. Two possibilities 
are therefore left open: a first order order phase transi- 
tion, or a second order phase transition with 0(4) critical 
indexes. 

Two completely different scenarios correspond to those 
two possibilities. Since, contrary to the first order case, 
second order phase transitions are unstable against the 
explicit breaking of the underlying symmetry, in the sec- 
ond case one would have a crossover instead of a real 
phase transition for small but non-zero quark masses. 
That would imply the possibility of going continuously 
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from a confined to a deconfined state of matter, in con- 
trast with the idea of confinement being an absolute 
property of strongly interacting matter at zero tempera- 
ture and of deconfinement being an order-disorder tran- 
sition associated to a change of symmetry [9| . A second 
consequence would be the presence of a crossover line, at 
finite mass, in the region of high temperature and small 
baryon chemical potential /is of the QCD phase diagram, 
thus implying a critical point fic| to connect with the first 
order line which is supposed to exist at low temperatures 
and large chemical potentials. No such point is expected 
to exist if the transition at (1b = is first order. No crit- 
ical point has been found up to now in ex per iments with 
heavy ions, but the question is still open [Til . [l2l. [l3j. 

The problem can in principle be solved numerically 
using lattice QCD simulations, by means of a finite 
size scaling (FSS) analysis leading eventually to a pre- 
cise determination of the critical indexes of the transi- 
tion. However technical difficulties are encountered in 
this program. First of all, huge computational resources 
are needed to make numerical simulations of the system 
which are close enough both to the chiral and to the ther- 
modynamical limit (i.e. with small enough quark masses 
and large enough spatial volumes). Second, a FSS analy- 
sis is made intricate by the fact that the system has two 
relevant scales: the correlation length £ of the order pa- 
rameter and the inverse quark mass l/m q . In particular 
it is possible to write the following FSS ansatz for the 
free energy density C/kT around the chiral critical point 

^^L-U(rL\l\a mq Ll-) . (1) 

C depends on two different scaling variables instead of 
one as is the case for simpler systems, like e.g. the 
quenched theory. L s is the spatial size, r is the reduced 
temperature t = (1 — T/T c ), v is the critical index of 
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the correlation length (£ ~ t~ v ) and yh is the magnetic 
critical index. From Eq. (JXJ) one can then derive the FSS 
for the specific heat 

C v - Co ~ L^ v tj> e (tL\I v \am q V>^) , (2) 

where Co stems from an additive renormalization [l4l |. 
and for the susceptibility x °f tl ic chiral order parameter 

Xm - xo - Vj v 4> x (rL x J v ,am q vf) . (3) 

A few groups have investigated the problem on the lattice 
with staggered pj, El El EI El Hi HI or Wilson Q 
fermions. The common procedure has been to use ap- 
proximate versions of the scaling laws and ((3]), usually 
assuming to be already in the infinite volume limit. No 
clear answer has been found, but there exists a general 
prejudice in favor of a second order chiral transition. 

In Ref. [l| we have approached the problem by use 
of staggered fermions and a novel strategy for the FSS 
analysis which, together with the availability of relevant 
resources of computer power, has allowed us to achieve 
some progress. We have decided to keep one of the scaling 
variable fixed, so as to reduce the FSS analysis, Eqs. @ 
and ([3]), again to one variable. In order to do that, after 
choosing a value for the critical index appropriate to 
a given universality class, we have performed a series of 
runs at variable spatial size L s and quark mass m q , keep- 
ing the quantity am q L v s h fixed. The choice for the index 
yh implies an a priori assumption about the critical be- 
havior, which can then be carefully checked without any 
approximation looking at the residual scaling. In partic- 
ular in Ref. |l[ we have chosen to test the 0(4) critical be- 
havior, or better 0(2), which is more appropriate for the 
case of staggered fermions at non zero lattice spacing [2(| : 
we have therefore fixed yh = 2.49, which happens to be 
the same both for 0(4) and for 0(2). Our results for 
the chiral susceptibility, for the specific heat and for the 
equation of state have then shown a clear inconsistency 
with both the 0(4) and 0(2) scaling hypothesis, thus 
giving clear evidence against the possibility that those 
critical behaviors can describe the QCD phase transition 
for N f = 2 in the chiral limit. In Ref. [J we did not 
perform the analogous scaling test assuming a first order. 
We do that in the present paper (see Sect. IIII|) . In Ref. [1] 
we have obtained some evidence in favor of a first order 
transition, by keeping the first scaling variable fixed. We 
then do the approximation of spatial size large compared 
to the inverse pion mass. With this approximation the 
scaling laws read 

C v - Co ~ am-^^Uc (tL 1 /^ , (4) 

Xm — Xo — am-^ v ^Ux {tL\/ v ) . (5) 

We have checked Eqs. (HJ)-© and we found disagreement 
with 0(4), 0(2) and agreement with a weak first order. 



A further result of our study was that a simple analy- 
sis of the dependence of the pseudocritical temperature 
on the quark mass, when correctly taking into account 
the dependence of the physical lattice scale on m q , does 
not allow to discriminate between a second order and a 
first order critical behavior. In Ref. [l[ also the magnetic 
equation of state was consistent with weak first order (see 
also Sect. Mil below). 

The results obtained in Ref. [l| must be considered 
as the starting point of an accurate study of the prob- 
lem. Indeed several questions have been left open, which 
deserve further analysis. First of all, the preliminary ev- 
idence in favor of a first order transition should be di- 
rectly confirmed by a series of runs in which am q L v s h is 
kept fixed according to first order critical behavior, i.e. 
yh = 3. Second, our evidence against a second order tran- 
sition in the universality class of 0(4) or 0(2) should be 
verified against all possible systematic effects which could 
have influenced our results, in particular: i) In Ref. [l[ we 
have used an non-exact R algorithm [23[, even if being 
quite conservative in the choice of the integration step for 
molecular dynamics; that could lead to systematic errors 
which, in principle, could influence the determination of 
the order of the phase transition, ii) In Ref. [l[ we have 
used a standard gauge and fermionic action, and a tem- 
poral extent L t — 4, corresponding to a lattice spacing 
a ~ 0.3 fm around the phase transition (T = l/(aL t )). 
Critical behavior is a typical infrared phenomenon. Nev- 
ertheless ultraviolet cut-off effects could in principle have 
some influence on it so that the continuum limit of our 
results should be checked by using a smaller lattice spac- 
ing, i.e. an improved action and/or a larger value of L t . 
Finally, if the chiral transition is really first order, how- 
ever weak, one should find signals of metastability in the 
physical observables when going to large enough volumes; 
no convincing signals were found in previous literature, 
nor in Ref. [1J, and the question should be clarified, by 
exploring larger volumes. 

Answering to all previous questions represents a diffi- 
cult and computationally demanding program, which we 
partially carry out in the present paper. In particular we 
address the question related to the use of a non-exact al- 
gorithm in Section [HI where some of the results obtained 
with the R algorithm in Ref. [l[ are checked by using an 
exact RHMC algorithm. In Section IIIII we directly test 
the first order hypothesis by using a new set of numeri- 
cal simulations performed by keeping am q L v s h fixed with 
yh = 3. Our conclusions and perspectives for the contin- 
uation of our program will be presented in Section HVl 



II. RESULTS OBTAINED WITH THE RHMC 
ALGORITHM 

Keeping under control the systematic errors introduced 
by a non-exact algorithm like the Hybrid-R is a very ex- 
pensive task in terms of computer power. In principle one 
should carry out multiple MC simulations with smaller 
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FIG. 1: Comparison between MC estimates obtained by the exact RHMC algorithm and the Hybrid-R for the average value 
(left column) and susceptibility (right) of the plaquette (top row) and chiral condensate (bottom) on lattices with L a = 16, 32 
and rriL — 0.01335. The average values for L s — 32 in the left figures have been shifted to the right (Af3 = 0.004) for sake of 
visual clarity. 



and smaller molecular dynamics integration steps for ev- 
ery value of the simulation parameters. This extrapola- 
tion is in practice unfeasible - especially if one aims at 
the investigation of the chiral limit. What can be done 
is to choose the integration step size as a function of 
the simulation parameters, in our case the lattice quark 
mass am q = , in such a way that discretization errors 
are negligible as compared to the statistical ones. The 
right functional form can be determined by a preliminary 
study on a representative subset of the parameters or us- 
ing known results present in the literature (see e.g. [2fj|). 
For standard Kogut-Susskind fermions the simple choice: 
St = rriL /4 is believed to lead to an accuracy of ~ 5% in 
the thermodynamic susceptibilities for masses as low as 
rriL ks 0.01. That was the choice we used in Ref. [l[ for 
all but the largest volume L s = 32 at the smallest mass 
rriL = 0.01335 where, for computational limitations, we 
took St — mi/2, which is expected to introduce errors of 
about 10%. 

Recently the RHMC algorithm [24[ has emerged as 
a convenient and exact algorithm for staggered fermion 



simulations. Not only it has no systematic errors to con- 
trol (i.e. no estrapolation is needed), but it also outper- 
forms the R algorithm for 2 flavors in terms of computer 
time. This makes the RHMC algorithm the ideal candi- 
date to put the data of Ref. [l[ to the test (and to produce 
the new data). 

The subset of simulations used for this comparison was 
the one with the smallest lattice quark mass, namely 
rriL = 0.01335 with L s = 16, 32, where discretization 
errors are expected to be more significant and where a 
larger St was used for the biggest lattice. We made three 
different simulations for each value of L s at couplings 
just below, just above and at the transition point. For 
each value of j3 we accumulated a statistics of order 3k 
trajectories for L s = 32 and of order 15k trajectories for 
L s = 16. 

The results are shown in Fig. Q] The two most sig- 
nificant quantities, the average plaquette and the chiral 
condensate, used in Ref. [l| are displayed together with 
their respective susceptibilities for both the lattices. 

A clear agreement at a ler level for all of these quan- 
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FIG. 2: Pseudocritical coupling curve in the (J3, tul) plane. 
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compared with Fig. 3 of Ref. [1( . 



tities is observed at the transition point. The other 
points agree within 2er except for the one correspond- 
ing to L s = 16 above the transition. The reason for this 
discrepancy can be traced back to the limited statistics 
used for the end tail in this simulation of Ref. [1], in a 
coupling region of small importance. In the scaling re- 
gion of interest the results obtained in Ref. [l[ are proved 
to be well within the expected errors. 

We confidently conclude that the use of the R algo- 
rithm in Ref. [l| did not introduce artefacts invalidating 
the finite-size scaling analysis. 



III. DIRECT TEST OF THE FIRST ORDER 
HYPOTHESIS 

Some hints of a first order transition were observed in 
Ref. 0, by use of Eqs. ®-@ using the dataset generated 
to directly check the 0(4), 0(2) universality classes. We 
now repeat the scaling analysis of Ref. [l[ assuming weak 
first order. We generate a dataset with a fixed value 
of am q L v s h , with — 3 as expected for a first order 
transition. The most convenient way to proceed is to 
use the already available and checked MC simulations at 
L s = 32 and tol = 0.01335. This enables a major saving 
in computer time. This also automatically fixes the value 
of ttilLs — 437.45. Three other sets of simulations were 
made to construct the dataset for the first order scaling 
test: one with L s = 16 and rriL = 0.1068, one with 
L s = 20 and tol = 0.054682 and one with L s = 24 and 
hil = 0.03164444. For each of them ten different values 
of [3 spanning the entire critical region were simulated 
with a total statistics of about 90k trajectories collected 
for each of the lattices. 

The pseudocritical couplings (3 C of the three new lat- 
tices are in excellent agreement with the pseudocritical 
curve determined in Ref.[l[ (see Fig. [5]). Given such an 



agreement, no modification is necessary to Sect. IVA of 
Ref.Q. 

As in Ref.yO] the background must be determined in 
order to check the scaling. For the case of the specific 
heat such background was observed not to depend on rriL 
and to be almost linear in in the mass region of rele- 
vance for our purposes. The three new datasets at mi = 
0.03164444, m L = 0.054682 and m L = 0.1068 nicely fit 
together with all other data. The background estimate of 
Ref.pl is only slightly modified by the new data (shifts 
of order 0.3cr): C (J3) = 0.417(51) - 0.0695(93)/3. For 
the chiral condensate susceptibility \m the background 
is determined by a best fit of the peaks of the curves to 
the function: x P eak(L s ) = xo + kL2 /v with j/v = 3, 
appropriate to a first order transition (see Eq. ([3])). 

The consistency check of the first order finite size scal- 
ing is shown in Fig. [3] (the analogous figures for 0(4) are 
Fig. 6 and 17 of Ref. p|). A reasonable scaling is observed 
for the specific heat Cy . As already stated in Ref. [l[ , Cy 
is independent of any prejudice on the symmetry and on 
the order parameter. Violations of the scaling Eq. ([3]) 
are observed for Xm a t larger values of the masses. In 
fact Eq. ([31) is expected to be valid for the susceptibility 
of the order parameter. At large masses chiral symmetry 
is badly broken and {ipifi) is possibly not a good order 
parameter. In our data, both of the present paper and 
of Ref. y|, Eq. §5§ seems to be violated for mj, > 0.05. 

An updated version of the scaling laws Eqs. (H)-© 
(Fig. 8, 9 and 18 of Ref. [J) with the new data is shown 
for completeness in Fig. 2] The scaling laws Eqs. (HI)-© 
correspond to taking the limit mi,L v a h 3> 1 with tL 



fixed in Eq.JT]). In the mass region considered a rea- 
sonable scaling is observed for first order for the spe- 
cific heat (similar considerations as above apply). For 
the chiral susceptibility good scaling is observed only for 
m L < 0.05. 

Both in Fig|3] and in Fig. 0] the background for the 
chiral susceptibility has been obtained by the fit to the 
curves with uil < 0.05. The fact that only in a neigh- 
boor hood of rriL = the scaling is expected is well known 
- but the actual mass range is not. The fitting procedure 
described allows us to identify a scaling window where 
the whole curve scaling can be verified. 

The scaling of the chiral condensate (magnetic equa- 
tion of state) can also be checked. The expectation for 
first order is: 



= /(t^l 1 ) 



(6) 



In Ref. [lj we found that if a subtraction is allowed on the 
left-hand side of Eq. © then a good scaling is found with 
the first order exponent. The subtraction was found to be 
numerically equal to the value of the condensate on the 
zero gauge field background (perturbative value). This, 
in the region explored, is equivalent to impose that all 
the curves coincide at the pseudocritical coupling. 

The updated figure for the first order scaling is shown 
in Fig. [5] (see Fig. 14 of Ref.[l|). The scaling is very good 
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within the statistical errors. The three new curves have 
been inserted and have been found to scale nicely. 

As a last check we looked for metastabilities in the MC 
histories of the new lattices. We found no metastable 
states. The history corresponding to the pseudo-critical 
coupling of the L s — 24 lattice is shown in Fig. [5] No 
clear double-peak structure is found. If the transition is 
first order, increasing the lattice size, one should eventu- 
ally see the metastabilities. At present, no clear sign of 
these states arc found in the numerical data. 



IV. CONCLUSIONS 

The determination of the order of the chiral transition 
for Nf = 2 QCD has turned out to be a very challeng- 
ing problem because of the huge computational resources 
required. 

In Ref. [l| we proposed a novel method to deal with 
the double scaling Eq. We applied such method 

to test the 0(4), 0(2) critical behavior - expected for a 
second order transition on the basis of theoretical specu- 
lations. The conclusion was that our numerical data was 
not compatible with 0(4), 0(2) and we found hints of a 
first order transition. 

In the present work we have completed our previous 
analysis in two different respects. 

In Sect. IIII1 we have performed a direct test of first 
order scaling, using new MC data with parameters, with 
m L L 3 s = 437.45, L s = 16, 20, 24, 32. These lattices show 
a good scaling for the specific heat and for the chiral sus- 
ceptibility at masses mj, < 0.05 (see Fig. [3]). At larger 
masses scaling of the chiral susceptibility is broken, pre- 
sumably because of strong breaking of the chiral symme- 
try. 



For the sake of completeness we have also updated the 
scaling pictures of Ref. [l[ corresponding to Eqs. ([I]),© 
(see Fig. |4j). 

The good scaling of the chiral condensate with first or- 
der pseudo-critical exponents, observed in Ref. is also 
in complete agreement with the new MC data produced 
for this work (see Fig. [5]) . 

In Sect, m we have checked some of the systematic er- 
rors present in the previous analysis of Ref. [l[ , by com- 
paring a representative subset of the old MC data gen- 
erated using the non-exact R-algorithm with a new one 
obtained with the exact RHMC algorithm. The direct 
comparison has shown no significant deviations between 
the two datasets at the critical coupling - thus validating 
the result presented in Ref. 

Taking all the evidence together, the first order scaling 
is clearly preferred over the second order 0(4) behavior. 
We can say that 0(4), 0(2) are excluded and first or- 
der is consistent with data, modulo possible effects due 
to the discretization. Again we think that ultraviolect 
effects should be irrelevant with respect to the large vol- 
ume behavior. However the use of finer lattices and new 
simulation algorithms to approach the chiral limit will 
possibly clarify this issue. 
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